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Abstract 

Background: Investigation of the nonlinear pattern dynamics of a reaction-diffusion system almost always requires 
numerical solution of the system's set of defining differential equations. Traditionally, this would be done by selecting 
an appropriate differential equation solver from a library of such solvers, then writing computer codes (in a 
programming language such as C or Matlab) to access the selected solver and display the integrated results as a 
function of space and time. This "code-based" approach is flexible and powerful, but requires a certain level of 
programming sophistication. A modern alternative is to use a graphical programming interface such as Simulink to 
construct a data-flow diagram by assembling and linking appropriate code blocks drawn from a library. The result is a 
visual representation of the inter-relationships between the state variables whose output can be made completely 
equivalent to the code-based solution. 

Results: As a tutorial introduction, we first demonstrate application of the Simulink data-flow technique to the 
classical van der Pol nonlinear oscillator, and compare Matlab and Simulink coding approaches to solving the van 
der Pol ordinary differential equations. We then show how to introduce space (in one and two dimensions) by solving 
numerically the partial differential equations for two different reaction-diffusion systems: the well-known Brusselator 
chemical reactor, and a continuum model for a two-dimensional sheet of human cortex whose neurons are linked by 
both chemical and electrical (diffusive) synapses. We compare the relative performances of the Matlab and 
Simulink implementations. 

Conclusions: The pattern simulations by Simulink are in good agreement with theoretical predictions. Compared 
with traditional coding approaches, the Simulink block-diagram paradigm reduces the time and programming 
burden required to implement a solution for reaction-diffusion systems of equations. Construction of the 
block-diagram does not require high-level programming skills, and the graphical interface lends itself to easy 
modification and use by non-experts. 

Keywords: Simulink modelling, Brusselator model, Cortical model, Turing-Hopf pattern 



Background 

Natural systems exhibit an amazing diversity of patterned 
structures in both living and nonliving systems [1], such 
as animal coats (e.g., zebra stripes, leopard spots and fish 
spirals), chemicals in a gel [2], laser light in a cavity [3], 
charges on the surface of a semiconductor [4], ecological 
balance between two species [5] and neuronal activations 
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in human cortex [6]. Generally, there exist two kinds of 
natural patterns [7] : 

• Thermodynamic equilibrium-based, such as crystal 
structures in inorganic chemistry and spontaneously 
emergent organic polymer patterns. The forming 
mechanisms of these patterns have been extensively 
studied and well-explained via thermodynamics and 
statistical physics: When a system's temperature 
decreases, its entropy and the Gibbs free energy 
accordingly become smaller (systems tend to move 
towards a thermodynamic equilibrium where the 
Gibbs free energy reaches a minimum). Thus the 



O 



© 2014 Wang et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
BiolVlGCl C6ITtr3l Commons Attribution License (http://creativecommons.Org/licenses/by/4.0), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication 
waiver (http://creativecommons.Org/publicdomain/zero/1 .0/) applies to the data made available in this article, unless otherwise 
stated. 



Wang etal. BMC Systems Biology 2014, 8:45 
http://www.biomedcentral.eom/1752-0509/8/45 



Page 2 of 21 



principle of minimum energy allows the system to 
form certain spatial structures. 
• Thermodynamics far-from-equilibrium, such as 
examples mentioned at the beginning of this section. 
These patterns emerge away from thermodynamic 
equilibrium, thus thermodynamic theory is not 
applicable. The understanding of these patterns may 
require application of kinetic theories. 

Pattern dynamics research focuses on universal pattern- 
forming mechanisms for the latter case. Away from ther- 
modynamic equilibrium, some experimentally observed 
2D spatial patterns have been reported: Rayleigh-Benard 
convection patterns [8]; reaction-diffusion Turing pat- 
terns [9]; Faraday- wave patterns [10]; vibratory granular 
patterns [11]; slime-mold spiral patterns [12] and Ca 2+ 
spiral waves in Xenopus laevis eggs [13]. Although these 
patterns are different in the spatiotemporal scales or 
detailed pattern-forming mechanisms, they have similar 
structures. 

More than half a century ago, British mathemati- 
cian Alan Turing strove to understand the spontaneous 
processes generating these structures. In his famous 
1952 paper [14], he proposed a reaction-diffusion model 
explaining potential mechanism for animal coats: At a 
certain stage of embryonic development, the reaction 
and diffusion between molecules, known as morphogens, 
and other reactors, lead to the breaking of symmetry of 
the homogeneous state. The morphogens spontaneously 
evolve to a non-uniform state, leading to the unique tex- 
tures seen on animal skin. The key pattern-forming idea 
in Turing's paper is the systems spontaneous symmetry- 
breaking mechanism, also known as a Turing bifurcation 
(spatial instability). 

However, for nearly 30 years Turing's paper drew 
little attention for two reasons: morphogens had not 
been found in biological systems; negative chemical 
concentrations are permitted in Turing's model, but 
are not accepted by chemists. From the late 1960s, a 
Brussels team led by Ilya Prigogine (winner of 1977 
Nobel prize in chemistry) was endeavouring to prove 
Turing's theory. Prigogine developed a reaction-diffusion 
model, Brusselator [15], to show the existence of Turing 
patterns that obey the rules of chemical kinetics. The 
Brusselator is a mathematical representation of the 
interaction between two morphogens, a reactor and an 
inhibitor, competitively reacting in time and spread- 
ing in space, which could give rise to a symmetry- 
breaking transition bifurcating from a homogeneous to 
patterned state, either stationary in a spatial (Turing) 
structure or in a temporal (Hopf) oscillation [15]. 
The Brusselator model further revealed the "secret" of 
Turing patterns: a coupling between nonlinear reac- 
tion kinetics and distinct diffusion rates such that the 



inhibitor should diffuse more rapidly than the activator 
[16]. 

This Brusselator proposed activator-inhibitor interac- 
tion is now a well known universal principle explaining 
regular pattern formation in chemical [17-19], ecological 
[5,20] and physical [21,22] systems, as well as biological 
systems [23,24]. 

As pattern-forming systems essentially consist of cou- 
pled differential equations, the simulated patterns are 
time- and space-dependent solutions of the differential 
equations. The MATLAB built-in fourth-order Runge- 
Kutta solver ode 4 5 and custom Euler methods are 
common ways to integrate differential systems. The 
implementation of these algorithms are well explained and 
demonstrated in many studies [25-28]. For example, in 
Yang's book [25], at the end of Part II Yang presents a piece 
of concise Matlab code for efficiently simulating sim- 
ple reaction-diffusion systems. With some modifications, 
Yang's programs can be used to simulate pattern forma- 
tion in a wide range of applications of nonlinear reaction- 
diffusion equations. Some of these examples are discussed 
in detail in Part III of Yang's book. Alternative to Mat- 
lab, there are other options for pattern simulations such 
as Meredys [29], implemented in the Java programming 
language, interpreting the NEUROML model description 
language [30]. Besides, there are other programming envi- 
ronments applicable to modelling pattern formations such 
as COMSOL [31] and Modelica [32]. 

Although it is efficient to solve differential equations 
in Matlab or other programming platforms, their code 
script-based pattern simulations require high-level pro- 
gramming skills to tune the model parameters in order 
to examine the pattern dynamics; this limits their uptake 
by non-programmers. A few attempts for graphic-based 
pattern simulations have been made in the last decades. 
For example, Ready [33] is an ideal program for getting 
started with reaction-diffusion systems. It runs fast (tak- 
ing advantage of multi-core architectures), is easy to use 
(graphic-based) and runs on multiple platforms. In addi- 
tion, there are other Java applets allowing easy pattern 
simulations such as Gray-Scott Simulator [34], showing 
the Gray-Scott pattern [35]: phenomena in grids of points 
that are connected "amorphously". This closely models 
the development of a biological system of living cells. 
Similarly, Lidbeck has created another Java application 
[36] with extensive presets and options, and even allows 
3-D of pattern simulations. However, these graphical 
user interface (GUI) -based softwares designed mainly for 
demonstrations of specific (pre-loaded) models. READY 
supports custom models but is written in READY-specific 
codes. So technically, Ready is still a coding-based plat- 
form with a GUI interface. In the community of pattern 
dynamicists, there may be a demand for a software plat- 
form circumventing the programming burden required to 
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implement numerical simulations of biologically-relevant 
pattern-forming systems. 

The aim of this paper is to introduce SlMULlNK mod- 
elling for simulating pattern dynamics. SlMULlNK, an add- 
on product to MATLAB, provides an interactive, graphical 
environment for simulating and analysing dynamic sys- 
tems. It enables modelling via a graphical programming 
language based on block diagrams. SlMULlNK has been 
used extensively in engineering field [37] such as control 
theory and digital signal processing for multidomain sim- 
ulation [38] as well as model-based design [39]. Besides, 
SlMULlNK users have extended its applications in other 
areas such as medical device prototyping [40], heat trans- 
fer modelling [41,42] and chemical reactions [43]. 

The present paper introduces the application of 
SlMULlNK to pattern simulation studies, and it is organ- 
ised as follows. We start, for pedagogical reasons, with 
a brief demonstration of the SlMULlNK in modelling dif- 
ferential equations. Then, we construct the well-known 
Brusselator model in SlMULlNK. By extending the Brus- 
selator modelling strategies, we construct a SlMULlNK- 
based human cortical model [44] (developed by the 
authors' team) that incorporates the pattern-forming 
essentials while retaining neurophysiological accuracy: 
the cortical model comprises interacting populations of 
excitatory and inhibitory neurons which communicate 
via chemical (neurotransmitter-controlled) and electri- 
cal (gap-junction) synapses. In the model, the interaction 
between chemical kinetics and electrical diffusion allow 
for the emergence of a comprehensive range of patterns, 
which may be of direct relevance to clinically observed 
brain dynamics such as epileptic seizure EEG spikes 
[6], deep-sleep slow-wave oscillations [45] and cognitive 
gamma-waves [46,47]. 

Finally, we examine the Brusselator and cortical model 
pattern dynamics. These simulation results are compared 
with the theoretical predictions. 

Methods 

Consider a generalised reaction-diffusion system for two 
reacting and diffusive species U and V of the form: 

d LI 

— =fu(U, V)+D U V 2 U 

8t (1) 

— =f v (U,V)+D v V 2 V 
at 

The diffusion constant Djjy [with units (length) 2 /time] 
is an important parameter indicative of the diffusion 
mobility. For a multi-component system, the higher the 
diffusivity, the faster the species diffuse into each other. 
Here,/^,y(Zi, V) is typically a nonlinear function of con- 
centrations U and V. 

Solving a pattern-forming system in the form of Eq. (1) 
requires interpreting the differential operators for time 



d/dt and space V 2 . In the following example, we first 
model the van der Pol oscillator in SlMULlNK to explain 
how we interpret the differential operator on time. 

Van der Pol oscillator 

The van der Pol oscillator was originally developed by the 
Dutch electrical engineer and physicist Balthasar van der 
Pol [48]. The van der Pol oscillator was the first mathe- 
matical model proposed for the heartbeat, and it has also 
been used to simulate brain waves [49,50]: 

^- / x(l-, 2 )-+, = 0 (2) 

We wish to solve this equation for the case fi = 1 with 
initial conditions x(0) = 2 and dx/dt = 0 at t = 0. 
The traditional way to solve a second-order differential 
equation is to convert to a pair of coupled first-order 
differential equations: 

x — y 

y = fi[l — x )y — x 

We would now integrate these equations with time 
using the MATLAB numerical integrator ode4 5. This 
helps to form the link with the integration in SlMULlNK. 

We code the first-order van der Pol equations into a 
MATLAB function a as follows: 



function dydt = vanderpoldemo ( t , y , Mu) 
%VANDERPOLDEMO 

%Defines the van der Pol equation for ODEDEMO. 

% Copyright 1984-2002 The MathWorks , Inc. 
% Revision: 1.2 Date: 2002/06/17 13:20:38 

dydt = [y(2) ; Mu* (1-y (1) "2) *y (2) -y (1) ] ; 



To solve Eq. (3), we specify the coefficient /x, the initial 
conditions and the time-span over which the integration 
is to proceed; then pass these values, along with the name 
of the van der Pol function, to the Runge-Kutta solver 
ode45: 



tspan = [0, 2 0] ; 
yO = [2 ; 0] ; Mu = 1; 
ode = @(t,y) vanderpoldemo (t , y, Mu) ; 
[t,y] = ode45(ode, tspan, yO) ; 

% Plot of the solution plot 
(t,y( : ,1) , t, y( : ,2) ) 
xlabel (' t' ) 
ylabel (' solution y' ) 

title ('van der Pol Equation, mu = 1') 



The calculated results are plotted in Figure 1. 

Alternatively, we may use the SlMULlNK construction of 
Eq. (2), as shown in Figure 2. 

At a first glance, the interface for SlMULlNK is com- 
pletely different from the code sheet. In SlMULlNK, all 
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van der Pol equation, ji = 1 




0 2 4 6 8 10 12 14 16 18 

Time [sec] 



Figure 1 Solution of the van der Pol equation, produced via Matlab code sheet. Program running time: 0.384 s in variable time-step. 
Simulation platform (same for all simulations in this paper): Matlab R201 3a, Mac OS X 1 0.9.1 , Xcode 5.0.2; CPU 2.4 GHz Intel Core i7, memory 8 GB 
1600 MHz DDR3. 



calculating elements are displayed by blocks. We select 
blocks from the SlMULlNK library, then connect them to 
build a model. 

The basic principle to model a differential equation in 
SlMULlNK is to find the input and output of an integrator. 
Since we have: 



d 2 x 



with the 



dx 
dt 



dt — x 



(4) 



then it follows that for a second-order differential 
equation, we need at least two integrators. As seen in 
Figure 2, we first place an integrator block (the left 
block labelled with j) to process the inner integration 

of Eq. (4): / ^?dt. The output of this integrator reads 
dx/dt, which is sent into the second integrator (the right 
block labelled with \). We assume the integrated x is 
known, thus being used to construct the input of the 



(5) 



left integrator block, which is equivalent to dt 
form: 

d x / o \ dx 

-—r = [1(1- X ) X 

dt 2 v 7 dt 

The product block (labelled with x) in Figure 2 com- 
bines (l — x 2 ) and dx/dt. The result is amplified by a gain 
(triangle block, valued /x), then passed through a function 
block where x is subtracted. Here, the RHS of Eq. (5) is 
constructed. 

Modelling a differential equation in SlMULlNK requires 
forming a closed loop, where the integrated variables are 
fed back into the system. Evolution proceeds until reach- 
ing the desired final time. The scope block shows the 
real-time output of the two integrators; the scope can be 
placed anywhere to monitor the response of a sub-system. 
The Outl and Out 2 terminals send outputs of two inte- 
grators to the Matlab workspace for further analysis. 
The results of this SlMULlNK model are exactly the same 
as shown in Figure 1. 



van der Pol Model 
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Figure 2 Si mu link built-in example for the van der Pol model called by the Matlab command vdp. 



Wang etal. BMC Systems Biology 201 4, 8:45 
http://www.biomedcentral.eom/1752-0509/8/45 



Page 5 of 21 



Both Matlab and Simulink allow fixed or self- 
adaptive (i.e., auto) time-steps for the Runge-Kutta 
solver b . Figure 3 shows that the discrepancy between 
Matlab and Simulink Runge-Kutta solvers in either 
fixed or auto time-step mode are sufficiently small 
(< 10 -10 ). Consequently, we can see that the accuracy 
of the model simulation does not depend on the mod- 
elling platform since MATLAB and SIMULINK share the 
same integration algorithm to solve differential equations. 
However, modelling in SIMULINK is more straightfor- 
ward and intuitive, and requires less programming skill 
than the Matlab code sheet. The original mathematical 
equations can be converted into SIMULINK by matching 
its pattern with Simulink blocks directly. Moreover, in 
Simulink, by simply adding more blocks, or replacing 
blocks, a new model is able to be built in a very short 
time. Simulink may be an ideal tool to efficiently per- 
form the simulations of a mathematical model. In the next 
section, we will extend the Simulink modelling method 
to describe a Brusselator system considering both its tem- 
poral and spatial evolutions on 1-D and 2-D Cartesian 
grids. 

Readers should be aware of the choice of an appropri- 
ate differential solver for a specific problem, depending on 
the stiffness of differential equations. Applying a wrong 
solver may lead to either unstable solution or exceptional 
computation time. However, it is practically difficult to 



identify the stiffness of a differential model, thus one 
should try at least two different solvers, and compare the 
results. If they concur, i.e. give the same solution, they are 
likely to be correct. As suggested by Matlab help file, 
it is worthwhile to try ode 4 5 first since it is the most 
widely used method. For pattern-forming systems, we can 
also compare the numerical solution with the theoreti- 
cal prediction to identify the applicability of the solver. 
For the demonstrated Brusselator and cortical models, 
ode 4 5 and ode 2 3 both work well and give rise to sim- 
ilar result; moreover, the numerical solutions match well 
with the theoretical predictions in emergent patterns (see 
Results and discussion section). So we choose ode4 5 
solver to integrate the differential-equation models in this 
paper. 

Brusselator model 

The Brusselator model describes the competition of 
two chemical species in a chemical reaction, and is the 
simplest reaction-diffusion system capable of generating 
complex spatial patterns. The competition between two 
reactors and the introduction of diffusion satisfy the key 
requirements for pattern formation [14]. The pattern 
dynamics of the Brusselator has been comprehensively 
examined in de Wit [15] and other researchers' work 
[51-54]. Here, our purpose is to introduce the 
Simulink modelling strategies. 



x 10 



(a) fixed time-step: 10 3 sec 




60 80 100 120 

Time [sec] 
(b) auto time-step 



200 




100 120 140 160 180 200 

Time [sec] 

Figure 3 Comparision for the van der Pol equation solved from Simulink and Matlab code sheet. Discrepancy over time for the solution of 
Eq. (2) is calculated from two modelling methods: Simulink and Matlab code sheet. Both methods use (a) fixed time-step 1 0 -3 s; or (b) auto 
time-step. Program running time for Simulink in fixed and auto time-steps are respective 1 .341 0 and 0.7404 s; for Matlab the corresponding time 
are 17.3952 and 9.4569 s. 
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Sim u li n k version of Brusselator model 

The simplest form of the model reads [55], 



—X = A - (B + 1)X + X 2 Y + D x V 2 X 
dt 

d 

— Y = BX-X 2 Y + D Y V 2 Y 
dt 



(6) 



where X and Y denote concentrations of activator and 
inhibitor respectively; Dx and Dy are diffusion constants; 
A is a constant and B is a parameter that can be varied to 
result in a range of different patterns. 

The LHS of Eq. (6) is a partial derivative on time since X 
and Y are functions of both time and space. At the RHS, 
the spatial derivative is represented by a Laplacian opera- 
tor V 2 . In the numerical simulation, the spatial dimension 
of the model is discretised into a grid by using the finite 
difference method. In the two-dimensional system the 
Laplacian with respect to the concentration field U in 
the node is calculated along the x and y directions 
simultaneously: 



A 2 



hi 



+ ■ 



h 2 



(7) 



where 



(8) 



The h Xf y demonstrators in Eq. (7) are the respective x 
and y grid spacings; they define the spatial resolution. 
Assuming h = h x = h y (i.e., a square grid), the dis- 
crete Laplacian operation in a one-dimensional Cartesian 
coordinates along the j-axis has the form: 



Ui,j+i — + Uij-i 

TP ; 



(9) 



for the two-dimensional case, we have 

Ui+y + Ui-y - W Uj + U Uj+ i + U Uj -i 



h 2 



(10) 



In Simulink, we initialise the Brusselator model as a 
column vector consisting of a 60 x 1 grid (spatial resolu- 
tion = 1 cm/grid-point) for the one-dimensional case; or 
as a 60 x 60 grid for the two-dimensional case. Grid edges 
are joined to give toroidal boundaries. 

The Laplacian operator V 2 D in Eq. (9) is implemented 
as a circular convolution of the 3x1 second-difference 
kernel L y 1D acting along the j-axis: 



'ID 



"ID 



h 2 



(ID 



The two-dimensional Laplacian operator V| D in Eq. (10) 
is built up from the sum of two orthogonal Lid operators: 



'2D ' 



"0 


1 0" 


1 


"0 


0 


0" 


1 


"0 


1 


0" 


0 


-2 0 


+ h 2 


1 


-2 


1 


'h 2 


1 


-4 


1 


0 


1 0 


0 


0 


0 


0 


1 


0 



(12) 



where we have again assumed a square grid so that h x = 

h y = h. 

In Simulink, the 1-D or 2-D Laplacian operator with 
toroidal boundaries is processed through two blocks: The 
"wrap-around" and "2-D CONV" (can process both 
1-D and 2-D convolutions) . The "wrap-around" block 
wraps the input matrix on both axes to allow a valid con- 
volution in the "2-D CONV" block against the Laplacian 
kernel L to return the cyclic convolution. We created 
a subsystem to compute the convolution, as shown in 
Figure 4. 

In Figure 4, the custom block labelled "wraparound" 
contains codes extracted from the convolve 2 ( ) 
function 0 . 



Out 



2-D CONV 



12 



V 2 X 



wraparound m 



MATLAB Function 



In2 



~\ In1 



In3 



Figure 4 Simulink modelling of the convolution with toroidal boundaries. The spatial derivative V 2 / approximates to the discrete 
convolution of /given by the kernal LX will be fed into the port lnl,the kernel L enters In2 and In3. 
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function y = wraparound (x, m) 








% Extend x so as to wrap around on both axe 




sufficient to allow a 


% "valid" convolution 


with m to return the 


cyclical 


convolution . 


% We assume mask origin near centre of mask 


for compatibility with 


% "same" option. 










[mx, nx] = size(x); 










[mm, nm] = size(m); 










if mm > mx | nm > nx 










error ('Mask does 


not fit inside array' ) 








end 










mo = f loor ( (1+mm) /2 ) ; 


no = f loor ( (1+nm) /2 ) ; 


% 


reflected mask origin 


ml = mo - 1 ; 


nl = no-1; 


% 


mask 


left /above origin 


mr = mm-mo; 


nr = nm-no; 


% 


mask 


right/below origin 


me = mx-ml+1; 


ne = nx-nl+1; 


% 


reflected margin in input 


mt = mx+ml ; 


nt = nx+nl ; 


% 


top of image in output 


my = mx+mm- 1 ; 


ny = nx+nm-1; 


% 


output size 


y = zeros (my, ny) ; 










y(mo:mt, no : nt ) = x; 


% central region 








if ml > 0 










y(l:ml, no : nt ) = x(me:mx, :); 






% top side 


if nl > 0 










y(l:ml, l:nl) = x(me:mx, ne:nx); 






% top left corner 


end 










if nr > 0 










y(l:ml, nt+l:ny) 


= x (me :mx, 1 :nr) ; 






% top right corner 


end 










end 










if mr > 0 










y(mt+l:my, no : nt ) 


= x ( 1 : mr , : ) ; 






% bottom side 


if nl > 0 










y(mt+l:my, 1: 


nl) = x(l:mr, ne:nx); 






% bottom left corner 


end 










if nr > 0 










y(mt+l:my, nt+l:ny) = x(l:mr, l:nr) 






% bottom right corner 


end 










end 










if nl > 0 










y(mo:mt, l:nl) = 


x ( : , ne : nx ) ; 






% left side 


end 










if nr > 0 










y(mo:mt, nt+l:ny) 


= x ( : , 1 :nr) ; 






% right side 


end 











The reason we introduce the custom block is that the 
SlMULlNK built-in 2-D CONV block provides only the 
"valid" (non-flux) boundary condition, and cannot handle 
periodic boundaries. 

Following the ideas of modelling the van der Pol oscilla- 
tor, we can easily convert Eq. (6) to SlMULlNK blocks, seen 
in Figure 5. 

Simulink versions of Waikato cortical model equations 

The authors' team at the University of Waikato (New 
Zealand) has developed a mean-field cortical model which 
encapsulates the essential neurophysiology of the cor- 
tex, while remaining computationally cost efficient [44]. 
The model envisions the cortex as a continuum in which 
pools of neurons are linked via chemical and electrical 
(gap-junction) synapses. 

Model background 

The Waikato cortical model has a rich history of develop- 
ment: The model is based on early work by Liley etal [56], 
with enhancements motivated by Rennie et al. [57] and 
Robinson et al. [58]; and has been recently extended to 
include electrical synapses (also referred as gap junctions) 



[44,47] supplementing neural communications via chem- 
ical synapses. 
The Waikato cortical modelling assumptions are: 

1. Cortical element is the "macrocolumn" containing 
-100,000 neurons. 

2. There are only two distinct kinds of cortical neurons: 
85% excitatory, and 15% inhibitory neurons. 
"Excitatory" and "inhibitory" are classified according 
to their effects on other neurons. 

3. There are three isotropic neuronal interactions: 
cortico-cortical (long-range from other 
macrocolumns), intra-cortical (short-range within 
macrocolumn) and connections from subcortical 
structures such as the thalamus and brain-stem. A 
macrocolumn with diameter ~1 mm and depth ~3 mm 
is sketched in Figure 6. We further assume that 
inhibitory connections via their short axons are 
locally restricted within a macrocolumn. In contrast, 
the axons from excitatory neurons are often longer 
and extensively branched, having length ranging 
from centimetres to several meters (e.g., in a 
giraffe's neck), allowing exclusively excitatory 
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Brusselator reaction-diffusion model (2D space) 



B X 



Pde1 



Initial state 



dX/dt 



Pde2 



Initial state 
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Figure 5 Si mu link construction of the Brusselator model. 
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cortico-cortical connections. Both excitatory and 
inhibitory connections are permitted in local neuron 
connections. 

4. Neurons exchange information via both chemical and 
electrical (gap-junction) synapses. Gap-junctions are 
more abundant within inhibitory populations, while 
being rare within excitatory neuronal populations. 

The authors' team first introduced gap-junctions into 
the cortical model in the paper Gap-junction mediate 
large-scale Turing structures in a mean-field cortex driven 
by subcortical noise [44] . The cortical model is expected 
to exhibit similar dynamics to a chemical reaction- 
diffusion system: The interaction between the excitatory 
and inhibitory neurons is analogous to to the competi- 
tion between the activator and inhibitor of the chemical 
reaction-diffusion model, e.g. the Brusselator. The gap- 
junction strength between inhibitory neurons also plays 
the same role as the diffusion terms in the Brusselator, 
which allows a spatial evolution of the patterns. Conse- 
quently, it is reasonable that the cortical model exhibits 
similar pattern dynamics for those seen in the chemical 
reaction-diffusion system. 



Model equations 

The neuron populations deliver spike flux 0^ from 
sources Q a at a distance following damped wave equations 
[58]: 

[(d/dt + vA ab f - v 2 V 2 ] </> ab = v 2 A 2 ab Q a (13) 

The subscript ab is read "a —> b" the connection from 
the presynaptic neuron type a to postsynaptic neuron type 
b. a, b can be either e (excitatory) or i (inhibitory). Q a is 
the spike-rate flux, which is a mapping from membrane 
voltage V a to population-averaged firing rates: 

Qmax 

® a = i + e -c(v a -e a )/a a ( 14 ) 

The total input flux arriving at the post-synapse is given 
as, 

M ab =NM+<A t + €1 (15) 

long-range local subcortical 

where superscripts a and f$ distinguish long-range and 
local chemical synaptic inputs; N* b and N^ b are the 
number of such input connections. The cortex is driven 
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|-< ~ 1 mm ►! 



Intra-cortical 




Figure 6 Schematic representation of the connective topology within a cortical macrocolumn. Only four of the ~1 00,000 neurons are 
shown. The shapes of neurons are determined based on their appearances under microscope: triangles are excitatory (pyramidal) neurons; circles 
are inhibitory neurons. 



by subcortical noise which enters the intra-cortex. Sub- 
cortical excitation 0^ is modelled as small-amplitude 
spatiotemporal Gaussian-distributed white noise super- 
imposing on a background excitatory "tone" (0^) whose 
constant level is set via a subcortical drive scale-factor 5: 

€b(f) = Web) + $eb(t)y/Iw%) (16) 

where % eb is the Gaussian-distributed white-noise sources 
being delta-correlated in time and space. Inclusion of 
white-noise stimulation is motivated by physiological evi- 
dence that the cortex requires a continuous background 
"wash" of input noise to function normally. 

The total input flux ® ab is the temporal convolution of 
the dendrite impulse response H (t) with the input flux 
M ab : 

<& ab = (dendrite response) ® (input flux) 

= H b (t)®M ab (t) 

= [ H b (t-t')M ab (t')dt' 
Jo 

where the dendrite dynamics are modelled by the alpha- 
function impulse response 

H b (t) = ylte-v* (18) 



Reducing y; is equivalent to increasing the time-to- 
peak {II yd for the hyperpolarising GAB A (gamma- 
aminobutyric acid) response, as indicated in Figure 7. 

By taking derivatives, Eq. (17) becomes 

(ji + y*) ®ab = YjM ab {t) (19) 

Therefore, the flux received by target neuronal popula- 
tions are: 

(it + Ye ) 2 = + ^ + ^ Y * (20a) 

Finally, those incoming spikes perturb the soma resting 
potentials: 

r b ^ = V r b eSt - V b + Pefeb^eb + Pifib^ib + Dl^ 2 V b 

chemical synapses electrical synapses 

(21) 

We write D bb as the diffusive-coupling strength between 
electrically adjoined e o e, i o i neuron pairs. To 
simplify the notation, we write (Di,D2) = (D ee ,Da). 

Symbol definitions for the cortical model are sum- 
marised in Table 1. 
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7z = 58.6 s" 1 

— 7i = 25 s" 1 



0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

Time [sec] 

Figure 7 Cortical IPSP (inhibitory post-synaptic potential) response (inhibitory alpha-function) curve. With delays in the inhibitory 
postsynaptic response (reduction of the inhibitory rate-constant yj), the occurrence of the IPSP peak will be delayed. Here, the position of the peak 
shifts from 0.01 7 to 0.04 s by reducing y\ from 58.6 to 25 s _1 . 



Sim u li n k versions of Waikato cortical model equations 

Let us first list the mathematical equations for the Waikato 
cortical model and examine their characteristics. 



• The cortico-cortical equation 

2 



• The intra-cortical equations 

'j t + *«* = [ N eb€„ + <Qe + €l] Ye 



Vv 2 



Kb = (vA e6 ) 2 Q, 



can be arranged by collecting temporal derivatives to 
the LHS: 

^r eh +2vA eb ^ h =v 2 V 2 4> a eh -v 2 A 2 4> a eh +(vA eb ) 2 Q e 



dt 2 



dt 1 



(22) 



have different RHS, but their LHS are in the same 
mathematical pattern: 

/a \ 2 a 2 a 9 



(23) 



Table 1 Symbol definitions for the cortical model 



Symbols 


Description 


Value 


Unit 




Inverse-length scale for e b axonal connection 


4 


cm -1 


V 


Axonal conduction speed 


140 


cm s _1 


synax 


Maximum firing rate 


30, 60 


s- 1 


&ej 


Sigmoid threshold voltage 


-58.5,-58.5 


mV 




Standard deviation for threshold 


3,5 


mV 


c 


Constant 


7T/V3 




N a eb 


Long-range e b axonal connectivity 


2000 






Local e b,i -> b axonal connectivity 


800, 600 




Ye 


Excitatory rate-constant 


170 


s- 1 


Yi 


Inhibitory rate-constant 


50 


s- 1 




e -> b tonic flux entering from subcortex 


300 


s- 1 


Tej 


Neuron time constant 


0.04, 0.04 


s 


e,i 


Reversal potential at dendrite 


0, -70 


mV 


v rest 


Neuron resting potential 


-64, -64 


mV 


AV S st 


Offset to resting potential 


1.5,0 


mV 


Pe 


Excitatory synaptic gain 


1.00 x 10~ 3 


mVs 


Pi 


Inhibitory synaptic gain 


-1.05 x 10~ 3 


mVs 


D 2 


/' ^> /' gap-junction diffusive coupling strength 


0-2.0 


cm 2 


Di 


e ^> e gap-junction diffusive coupling strength 


D 2 /100 


cm 2 
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We can move the term y 2 0 to the RHS of the 
intra-cortical equations, then the LHS of the 
intra-cortical equations have the expression: 

+ (24) 

which is similar to the LHS of the cortico-cortical 
equation. 
• The soma equation 

T&— - = V r b est -V b +(pe^eb^eb+Pifib^ib)+D bb V 2 V b 

01 

can be re-arranged as 



dt 



T b 



[Vl eSt ~Vb + (Pefeb®eb+Pitib®ib) +D bb V 2 V b ] 

(25) 



Following the ideas of SlMULlNK modelling in van der 
Pol oscillator, we need two integrator blocks for Eqs. (22) 
and (24), and two convolution processing for Eqs. (22) 
and (25). 

The strategy for modelling a large system is to focus 
on its subsystems first, then connect them together. The 
Waikato cortical model has three major parts: cortico- 
cortical, intra-cortical and soma equations. Figure 8 shows 
how neuronal fluxes are transferred from one to another: 
cortico-cortical flux 0^ is delivered to the long-range 
targets O ee and intra-cortical flux $> ee and 
0/ e and ®u merge into the soma equations. The out- 
put of the soma V e is connected to source neurons to 
form the closed loop through the excitatory sigmoid 
function: 

fjmax 

Qe = 



1 -f e -C(V e -O e )/<T e 



In following sections, we detail the SlMULlNK imple- 
mentation of the three subsystems (drawn as three blocks 
in Figure 8) of the cortical model. 



Cortico-cortical flux The SlMULlNK based cortico- 
cortical block (see Figure 9) is converted from Eq. (22). 
The flux-source Q e is a mapping from the excitatory soma 
voltage sent via the Goto block, to the firing-rate received 
via the From block. After two integrations, signals will 
be passed to the excitatory synapses via port-1 (upper 
right corner) and inhibitory synapses via port-2 in the 
intra-cortex. 



Intra-cortical flux In SlMULlNK modelling, we divide 
Eq. (16) into two parts: the constant (i.e. dc-level) excita- 
tory background s(<p^} and the one-off kick 
As illustrated in Figure 10, we use a Clock block to count 
the iteration step. Once the counter is above one, the 
"switch" will turn off the kick, allowing only the constant 
excitation to enter the intra-cortex (removing the Clock 
block would allow on-going noise stimulus from the sub- 
cortex). The 2-D spatial white noise are generated by the 
Band-Limited White Noise block. 

The intra-cortical model describes how post-synaptic 
fluxes evolve over time. In Figure 11, the local fluxes (input 
via the From block) along with the long-range fluxes 
(from input- 1 at the left, labelled as 0") and subcorti- 
cal drive are summed, then filtered at the post-synaptic 
dendrite, thus forming the post-synaptic fluxes O ee at the 
output port-1 (upper right corner). The O ee and O e ; flux 
models have symmetric structure. 

We assume that the cortico-cortical fibres are exclu- 
sively excitatory, thus there are no long-range inhibitory 
fluxes entering into the soma. Figure 12 shows that 
the local inhibitory fluxes <$>; e come from local source 
Qi only. The 4>/ e and models also have symmetric 
structure. 



Soma voltage Figure 13 presents the soma model of 
the excitatory neuronal group. The short-range fluxes 
are accumulated at the soma, from ports 1 and 2. The 



sigmoid mapping 



cortico-cortical 



(...)<l?ee = [iVe a e C + ---]7e 
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Figure 8 Flux flows for the Waikato cortical model between its cortico-cortical, intra-cortical and soma equation subsystems. 
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Figure 9 Si m u li n K-based cortico-cortical wave-equation. 



soma voltage V e is converted to firing-rate Q e locally 
in this sub-model (block labelled with Q e sigmoid), 
then fed back into the cortico-cortical and intra-cortical 
models. 

Finally, we connect all subsystems to form the com- 
pleted Waikato cortical model, as illustrated in Figure 14. 
It follows the flux flow-chart of Figure 8, with the 
detailed Simulink block connections shown in Figure 15. 
We argue that such model-based-design is an advan- 
tage for representing differential equations in SIMULINK. 



Although Simulink is useful for rapid prototyping, the 
Simulink implementation of the cortical model runs 
slower than our pre-coded Euler integration d since it 
is time-consuming to interpret the MATLAB function 
wrapround (see Figure 4) in SIMULINK. A 60 x 60 grid 
(side length 20 cm) 5-s cortical simulation takes ~10 s via 
Matlab Euler integration (fixed time-step 0.8 xlO -3 s), 
while ~40 s via Simulink (auto time-step and in accel- 
erator mode). Thus, it is recommended to avoid using 
Matlab functions in Simulink unless necessary. 
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Figure 10 Si mulin K-based sub-cortical flux. 
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Figure 11 SiMULiNK-basede e post-synaptic flux & ee for the intra-cortex. 
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Results and discussion 

Brusselator pattern simulations in 1-D space 

Before simulation, we first predict the pattern dynamics 
of the Brusselator model. Following the work by Steyn- 
Ross etal. [44], we applied a linear stability analysis (LSA) 



[59] by examining the sign of the real and imaginary parts 
of the dominant eigenvalue (wavenumber q dependent) 
o q = a(q) + ico(q) derived from the Jacobian matrix at 
a steady-state. LSA states that the stability of the steady- 
state: a(q) > 0 leads to unstable steady-state; a(q) < 0 
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Figure 12 SiMULiNK-basedi e post-synaptic flux <I>/ e for the intra-cortex. 
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Figure 13 SiMULiNK-based excitatory soma equation. 
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leads to stable steady- state. The scaled imaginary part 
Q)(q)/2n predicts the oscillation frequency of the pattern 
at the wavenumber q. Considering different combinations 
of the sign of a and co, there are four main classes of insta- 
bility, summarised in Table 2. Following these rules, the 
LSA shown in Figure 16(a) and (b) predict respectively a 
Hopf and Turing instability. The simulation in the upper 



panel shows a homogeneous oscillations, in the bottom 
panel exhibits a frozen spatial structure. Both simulations 
are in good agreement with the LSA predictions. 

Brusselator pattern simulations in 2-D space 

Brusselator simulation in a 2-D space has more vari- 
eties than in a 1-D space. These 2-D patterns will have 
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Figure 14 Overview of the Si mu link implementation of the Waikato cortical equations. 
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Figure 15 Detailed connection diagram for the SiMULiNK-based Waikato cortical equations. Solid arrow: direct connection; dashed arrow: 
Goto-From connection. 



unique spatial Turing structures that we cannot infer 
from the LSA. To precisely predict the Turing pattern 
dynamics, we derived the amplitude equations [54] for 
the Brusselator model at the onset of a Turing instabil- 
ity. The amplitude equation describes a reduced form 
of a reaction-diffusion system yet still retains its essen- 
tial dynamical features. By approximating the analytic 
solution, the amplitude equation allows precise predic- 
tions of the pattern dynamics when the system is near 
a bifurcation point. In simple words, the amplitude 
equations may offer a guidance for model parameter 
settings corresponding to specific patterns, e.g., Turing 
pattern textures. 

An application of the multiple-scale expansion [60] 
on the Brusselator model yields the following amplitude 
equation for the Turing mode [61]: 

^-Z^nZ^vZ^Zt-g \Z 1 \ 2 Z 1 -h(\Z 2 \ 2 + \Z 3 \ 2 ) Zi 

(26) 



where 



fi = (B-B c )/B Cf 

_ 38^^+5(A^) 2 -8-8(A^) 3 
g — 9A 3 r)( l+Ar)) 
TJ = ^Dx/Dy 



2 (I^Ar) , 2 
V — A \1+Ar] ) + A^ f 

h 5^+7(^) 2 -3-3(^) 3 
A 3 ri(l+Ari) 



in which B c is the critical value (see [54] for its set- 
ting) for Turing condition. Here, Z(i,2,3) describes slow 
modulations of the Turing pattern, in short we call it 
Turing amplitude. The equations for Z2 and Z3 can be 
obtained from Eq. (26) by simple permutation of the 
indices. Indexed subscripts indicate that three amplitudes 
correspond to their respective wavevectors qi, #2 and #3 
oriented at 120° angular separations (see Figure 17). Ide- 
ally, the resonance of the three wavevectors result to a 
hexagonal or stripes modes. Practically, these modes may 
mix to form distorted patterns. 

Following the work by Verdasca et al [61], we applied 
linear stability analysis to the steady-state solutions of 
Eq. (26), leading to the pattern prediction diagram, shown 



Table 2 The onsets of four main classes of instability 



Pattern stability 



Bifurcation 


Critical wavenumber 


Eigenvalue [a q = a + ico) 


Spatially 


Temporally 


Turing 


q^O 


a q =0 


Unstable 


Stable 


Hopf 


q = 0 


a = 0, co ^ 0 


Stable 
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a = 0,co^0atq = 0 


Unstable 


Unstable 


Wave 
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Figure 16 Brusselator simulation in 1-D space. Left column: dispersion curve of real (a) and imaginary (&>) parts of dominnat eigenvalues 
predicting the emergent pattens for two sets of parameters: (a) A = 2.5, B = 9, D x = 7, D Y = 1 0; (b) A = 2, B = 4.8, D x = 2, D Y = 1 0. H: Hopf 
mode with a > 0, co > 0 at wavenumber q = 0; T: Turing mode with a > 0, co = 0 at q ^ 0; DT: damped Turing with a < 0 at q ^= 0; DH: damped 
Hopf with a < 0 at q = 0. Right column: one-dimensional Brusselator model of length 60 cm with periodic boundary condition evolves in time 
running rightwards during 30 s. Colour indicates the local concentration of the reactant: [red] high concentration, [blue] low concentration. 



in Figure 18. The diagram comprises of three stability 
curves revealing the stability of specific Turing modes: 
a honeycomb-like hexagonal structure, stripes, and reen- 
trant honeycomb. Verdasca et at. denote the hexagonal 
mode as and its reentrant form as Hq (see [61] for 




Figure 17 Wave vectors of a hexagonal pattern. Superposition of 
three wave vectors at an angle of 1 20 degree with each other to form 
a hexagonal pattern. 



Verdasca etal/s explanation on the subscripts n and 0). In 
this study, \i is a bifurcation control parameter, the value 
of which leads to a stable (solid curve) or unstable (dashed 
curve) mode. 

For the 2-D spatiotemporal simulation of the Brusse- 
lator model, Simulink was set to 50-s simulation time 
in auto time-step mode, allowing the pattern-forming 
system to organise itself sufficiently. Figure 18 predicts 
the stabilities of stripes (red), Ho (blue) and (black) 
modes for the Turing instability of the Brusselator model. 
From the range of solid curves, we have the summary 
of parametric space where a specific mode is stable: The 
stripes mode is stable when fi~ < fi < fif; the hexagonal 
mode is stable (H^ and Ho) when fi < fi^ or fi > /XjJ". H^ 
and Ho interact at /x p where they exchange mode stability, 
that is, H^ will transit to Ho when \i crosses /x p from its 
LHS to the RHS. 

To verify the predictions from the bifurcation diagram, 
we select five different values of fi then examine the 
simulated patterns: 

(a) fi = 0.0495 falling into a range where only H^ is 
stable; 

(b) fi = 0.1100 falling into a range where stripes and H^ 
modes coexist; 

(c) /jL = 0.3994 where only the stripes structure is stable; 

(d) [A = 0.7000 again falling into a bistable range where 
stripes and Ho coexist; 

(e) fi = 1.4802 where only Ho is stable. 
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Figure 18 Turing mode stability of the Brusselator model in 2-D space. Each coloured stability curve represents specific mode: red = Stripes, 
blue = H 0 , black = . Solid and dashed curves correspond to stable and unstable modes respectively, according to mode stability analysis. Five 
representative /x values are selected for comparison of theoretical predictions for mode stability against practical simulations (shown as subplots). 
Colour of the pattern indicates the local concentration of the reactant: [red] high concentration, [blue] low concentration. Model parameters: 
A = 5, D X = 5,D Y = 40. 



In Figure 18, we see good agreement between 
SlMULlNK simulated patterns and theoretical predictions. 
Clear H^, stripes and Ho structures are observed at (a), 
(c) and (e) cases respectively, (b) and (d) show mixed 
states between forward and backward stable structures. 
Consequently, we can conclude that by increasing /x, the 
distance to the critical point, positively, Brusselator forms 
sequentially from H^, to stripes, to Ho. 

In summary, to construct the Brusselator model in 
SlMULlNK, we first place an integrator block to repre- 
sent time derivative. The temporal integrators output will 
be fed back into the system to engage with the systems 
evolution, then form the input of this integrator block, 
closing the loop. The model parameters can be adjusted 
by tuning the settings of the blocks A and B as well as 
two gains (labelled Dx and Dy respectively). The real-time 
spatiotemporal evolution of X and Y are monitored via the 
Matrix viewer block. The simout block delivers the 
solution of Eq. (6) to the MATLAB workspace for future 
analysis. The solution is a three-dimensional matrix with 
the third dimension the same length as the time span. 

In the next section, we will implement the SlMULlNK- 
based cortical model and examine its clinically-relevant 
pattern dynamics. 

Cortical model stability and simulations 

The work by Steyn-Ross et al. [44] suggests that the 
strong gap-junction diffusivity (large D2 value in the 
model equation (21)) provides a natural mechanism 
for Turing bifurcation that leads to the spontaneous 
formation of Turing labyrinth patterns of high and low 



neural activity that spread over the whole cortex, allowing 
multiple, spatially separated cortical regions to become 
activated simultaneously. Figure 19 illustrates the emer- 
gence of such cortical Turing patterns provoked by a 
strong inhibitory diffusion D2. It is possible that the 
Turing spatial synchrony explains the cognition "binding" 
phenomenon, which is, widely separated neural popu- 
lations that are anatomically unconnected are in very 
similar states of activity, thereby becoming functionally 
connected and giving rise to coherent percepts and 
actions. 

At the vicinity of a Turing instability, a weak Hopf 
instability can be induced in parallel by prolonging the 
timing of delivery of inhibition at chemical synapses. This 
permits slow Hopf oscillations with the spatial structure 
maintained. Specifically, a reduction of the inhibitory rate- 
constant yi in Eq. (20b) below a critical value ~30.94 s _1 
is sufficient to produce a complex dominant eigenvalue 
at zero wavenumber whose real part is positive; thus 
suggesting a global Hopf oscillation. 

In Figure 20(d), setting y { = 29.45 s" 1 predicts a ~0.95 
Hz Hopf oscillation. Independently, a Turing instability is 
boosted with moderately strong inhibitory diffusion D2 = 
1 cm 2 above its critical value 0.9066 cm 2 . Cortical Turing- 
Hopf interactions lead to beating patterns revealed in the 
time series recorded in Figure 20(b) for a single pixel on 
the cortical grid. The Fourier spectrum shows two fre- 
quency components whose difference matches with the 
ultra-slow envelope frequency, which is likely to be the 
weakly-damped resonance at ~0.152 Hz. 

Steyn-Ross et al. [46] posit that the interacting low- 
frequency Hopf and Turing instabilities may form the 
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(a) Time = 0 sec (b) Time = 0.9 sec (c) Time = 1 sec (d) Time = 3 sec 




Figure 19 Simulated cortical Turing patterns. Four snapshots (a-d) taken from a 3-s simulation of the cortical Turing pattern of the excitatory 
firing-rate Q e commencing from a homogeneous equilibrium. The top panel is the 3-D Q e plot, and the bottom panel is the same information 
viewed on a bird's-eye image. The gap-junction diffusity Dj = 1 cm 2 . The cortical model is initialised as a 1 00 x 1 00 grid with 20-cm side-length. 
Simulink is set to auto time-step mode. Simulation running time ~30 s. 



substrate for the cognitive state, namely, the "default" 
background state for the non-cognitive brain. These slow 
patterned oscillations may relate to very slow (<0.1 Hz) 
fluctuations in BOLD (blood-oxygen-level dependent) 
signals detected using fMRI (functional magnetic res- 
onance imaging) of relaxed, non-tasked human brains 
[62,63]. Steyn-Ross etal. [47] also predict that this default 



state will be suppressed with elevated levels of subcortical 
drive during goal-directed tasks [64-68]. 

The effect of embedding Matlab functions in 
Simulink running efficiency 

Although Simulink has an intuitive programming logic 
and comparable accuracy to Matlab, it sometimes runs 



cortical interacting Turing and Hcpf instabilities 



eigenvalue dispersion curve of the cortical model 




1.5 

Frequency [Hz] 

Figure 20 Beating-wave patterns of the cortical model. With strong gap-junction diffusion D 2 and carefully chosen inhibitory rate-constant y h 
eigenvalue dispersion curve (d) predicts a mixed pattern of Turing and Hopf instabilities, or and co are the real and imaginary part of the dominant 
eigenvalue respectively. Through a 400-s Simulink simulation, the Fourier spectrum (c) indicates a 0.1 5-Hz ultra-slow oscillation of the beating 
pattern (b) zoomed from (a) 0 e time evolution of the point at position (1 , 30) shown in (e) 25-x25-cm grid (1 00 x 1 00 grid-points) 3-D plot. 
Simulation running time ~75 min. (Figure modified from [46]). 
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much slower than Matlab, e.g., in the demonstrated 
Brusselator and cortical model simulations. The reason 
is that we embed MATLAB functions wraparound in 
the model to expand the SlMULlNK capability. Once 
a Matlab function block is present, the Matlab 
interpreter is called at each time-step. This drastically 
reduces the simulation speed. So, one should use the 
built-in blocks whenever possible. Without using Mat- 
lab function blocks, SlMULlNK shows a higher per- 
formance than Matlab, e.g., see the description of 
Figure 3. MathWorks Support Team also presented com- 
prehensive guidance to speed up the SlMULlNK simula- 
tion, which are available at http://www.mathworks.com/ 
matlabcentral/answers/94052. In the further optimisation 
of our SlMULlNK model, we will consider replacing MAT- 
LAB function with the MEX S -function, which may help 
to accelerate the simulation in the merit of its direct 
communication with the SlMULlNK engine (avoid the 
time consuming compile-link-execute cycle). The build 
of MEX S -function requires higher level programming 
skills, so we only address a more accessible method 
(embedding MATLAB function) for non-experts in this 
paper. 

Conclusions 

The strategy for modelling differential equation models 
in SlMULlNK is addressed in this paper. To construct 
a system of differential equations in SlMULlNK, we can 
directly convert the mathematical terms and operators 
to the SlMULlNK graphical block diagrams. The key idea 
for programming differential systems in SlMULlNK is to 
form a closed loop, such that the solution of the system 
can evolve in this loop. The accuracy and reliability of 
the SlMULlNK modelling method has been examined via 
comparing a van der Pol oscillator represented in Mat- 
lab code-script and SlMULlNK block-diagram, showing 
the SlMULlNK model has comparable performance with 
the code-script version. 

Using a well-known Brusselator model, we demon- 
strated two main SlMULlNK modelling strategies for a 
reaction-diffusion system: interpretation of 2-D convolu- 
tion with the periodic boundary condition; hybrid pro- 
gramming in SlMULlNK with MATLAB functions. The 
pattern simulations of the Brusselator in SlMULlNK are 
in good agreement with the predictions via bifurcation 
theories. 

Following the pattern-forming theories, we introduced 
a cortical model with competitive neuronal interactions 
and diffusions. Unlike the simple Brusselator, the cor- 
tical model is a complicated system comprised of dis- 
tinct cortical connections. Here, we showed how to build 
these subsystems in SlMULlNK. Finally, we connected 
all subsystems to form the completed Waikato corti- 
cal model. The simulations of the cortical model exhibit 



Turing and mixed Turing-Hopf patterns, the clinical rel- 
evance of which is briefly discussed. As the main aim 
of this paper is to introduce SlMULlNK modelling, read- 
ers are referred to Steyn-Ross et al.'s recent publications 
[6,44-47] for comprehensive investigations of the cortical 
model. 

Endnotes 

a vanderpoldemo is a MATLAB pre-coded function 
b fixed time-step ODE solvers are not built into 

MATLAB, but they can be acquired from a release by the 

MathWorks Support Team: 

http://www.mathworks.com/matlabcentral/answers/ 
uploaded_files/5693/ODE_Solvers.zip 

c The two-dimensional circular convolution algorithm 
was written by David Young, Department of Informatics, 
University of Sussex, UK. His convolve 2 ( ) code can 
be downloaded from MathWorks File Exchange: 
http://www.mathworks.com/matlabcentral/fileexchange/ 
226 1 9- fast- 2- d- convolution 

d MATL AB simulation codes were written by Alistair 
Steyn-Ross. The complete codes, plus README files and 
movies of cortical dynamics, are available from the web 
site: http://www2.phys.waikato.ac.nz/asr/ 
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